Journal of Computational Biology

Constructing microbial association networks is a common strategy for exploring relationships among taxa in microbiome studies. Although marginal correlation methods are easy to implement and allow formal inference, they can produce spurious edges driven by indirect associations through other taxa. Conditional graphical-modeling methods aim to recover direct associations, but many rely on Gaussian or linear assumptions and often provide limited uncertainty quantification. We propose a conditional, nonparametric approach based on the scaled expected conditional covariance (SEcov). SEcov measures population-level conditional association by residualizing each taxon with respect to the remaining taxa and scaling the resulting expected conditional covariance. The resulting estimator can incorporate flexible machine-learning methods for conditional-mean estimation and admits asymptotic normal inference, enabling p-values and confidence intervals for taxon-pair associations. We demonstrate through simulation studies that our proposed approach improves network recovery relative to other methods, and we illustrate the new method via construction of a co-occurrence network for the vaginal microbiome during pregnancy. IMPORTANCEHigh-throughput sequencing has made it possible to characterize microbial communities at large scale, and network analysis is widely used to summarize relationships among taxa. However, networks based on marginal correlations may include indirect associations, whereas many conditional graphical models rely on assumptions that may be difficult to justify for sparse, zero-inflated, compositional microbiome data. SEcov offers a practical alternative by estimating conditional associations nonparametrically and attaching inferential uncertainty to individual edges. This allows investigators to construct microbiome networks using statistically interpretable evidence for taxon-pair associations, rather than relying solely on arbitrary correlation cutoffs or regularization tuning parameters.

Motivation: A continued decrease in sequencing costs has facilitated the exponential increase in available sequencing data, with public databases like the European Nucleotide Archive (ENA) and Sequence Read Archive (SRA) reaching well in the order of petabases. This has been the incentive to develop more scalable tools for common bioinformatics tasks. One such task is the approximate searching of short sequence patterns like genes or reads in reference data sets. In recent years, a variety of indexing data structures have been proposed for searching large sequencing databases. The state-of-the-art index, the Hierarchical Interleaved Bloom Filter (HIBF) was first-in-class to index one million samples. To be useful for expanding repositories, it must be extended to support dynamic updates. Results: In this paper, we introduce a scalable and updatable sequence-search index by extending the HIBF with partial rebuilding to support efficient updates. We demonstrate the Dynamic HIBF's capacity for large-scale data by iteratively creating an index from over 100 TB of compressed reads across more than 39,000 full human RNA-Seq samples, updated in consecutive batches of 100. To benchmark against state-of-the-art tools, we evaluated incremental performance on a subset of 5,000 samples sub-sampled to 1% of their original read depth. In this comparative setting, the dynamic HIBF completed the sequential insertion of all 5,000 samples within 5 hours--24 to 65 times faster than competing methods and twice as fast as the static HIBF.

We describe LRSPAT (low-rank spatial toolkit), a fast and memory-efficient framework for approximating measures of spatial association for high-dimensional data. While LRSPAT can be applied to any multivariate spatial dataset, development was motivated by the computational challenge of identifying spatially variable genes in high-resolution spatial transcriptomics (ST) data generated by technologies such as 10x Visium HD, Xenium and Atera. LRSPAT leverages a truncated SVD of the expression data and a thresholded spatial weights matrix to perform reduced-rank reconstruction of spatial statistics in the quadratic form family, including global and local versions of Moran's I, Geary's C, and Getis-Ord G. A regularization approach is leveraged to account for the inflated null distribution of spatial statistics computed on latent variables. By performing key operations on the low-dimensional embeddings, LRSPAT is orders of magnitude faster than standard implementations with significantly lower memory requirements. Because the low-rank approach denoises and desparsifies ST data, LRSPAT is also more accurate than standard techniques at identifying genes with true spatial expression patterns. The dramatic improvements in execution time and memory consumption enable the genome-wide analysis of spatially variable genes (SVGs) and exploration of the full range of hyperparameters including spatial scale, distance metric, and embedding rank. This preprint outlines the background and mathematical details of the approach with limited preliminary results and a short conclusion.

Standard methods for ancestral sequence reconstruction (ASR) rely on substitution models for the residues in a biological sequence and assume independent evolution across these sites, ignoring the epistatic interactions that shape molecular evolution. In contrast, deep learning models like variational autoencoders (VAEs) can learn low-dimensional representations ("embeddings") of sequences in a protein family that may implicitly handle these dependencies, raising the possibility of performing more accurate ASR by interpolating between extant sequence embeddings within the VAE's latent space. In this study, we test this hypothesis by developing and evaluating a VAE-based ASR pipeline. Benchmarking this approach against established likelihood-based and parsimony methods using various simulations of protein evolution, including scenarios with and without epistasis, we find that the VAE-based approach is consistently and significantly outperformed by standard methods, even in epistatic regimes where it was hypothesized to have an advantage. We further show that this failure is not due to a lack of phylogenetic structure in the latent space, which does contain evolutionary signal. Rather, the primary limitation is the information loss inherent to the autoencoding process: the VAE's decoder cannot generate sequences with sufficient fidelity for the precise demands of ASR.

THAP9 is a transposable element-derived gene which encodes a protein that is homologous to the active Drosophila P-element transposase (DmTNP). Both THAP9 and DmTNP possess a C-terminal domain (CTD) which is functionally uncharacterized. Sequence and structural analysis suggest that the THAP9-CTD has a novel fold which is only found in THAP9 homologs. To explore the evolutionary history and characteristics of this novel domain, exhaustive phylogenetic analysis (using MSA, structure prediction, MSTA-based clustering) was performed. THAP9-CTD homologs were more widely distributed throughout the animal kingdom in comparison to DmTNP-CTD homologs which were restricted to arthropods. Moreover, the THAP9-CTD homologs were more conserved, especially among mammals and birds and their average length increased in a class-specific manner. Comparison with the DmTNP-CTD homologs demonstrates that although their respective CTDs may have evolved independently, they both surprisingly share similar secondary structure elements consisting of three conserved helical regions made of hydrophobic residues that are predicted to make up a conserved core. The role of the respective CTDs were further investigated by creating truncation mutants lacking the CTD. Interestingly both THAP9 and DmTNP truncation mutants are still capable of DNA excision and integration suggesting that their respective CTDs are not essential for DNA transposition. Moreover, CTD truncation favours DNA integration in THAP9: this suggests that CTD acquisition during evolution may have led to THAP9 domestication as observed in other transposable element-derived genes like Rag1 and piggybac, which have similar terminal regulatory domains.

Bioinformatics software and databases are essential components of modern life science research, yet their mentions in the scientific literature are often inconsistent and difficult to systematically identify at scale. The lack of a comprehensive and up-to-date catalog of bioinformatics resources hinders efforts toward automated biomedical knowledge extraction and streamlined data analysis. Here we present SNAIL, a hybrid named entity recognition framework designed to automatically identify bioinformatics software and database (SW/DB) names from biomedical texts. SNAIL integrates complementary lexical and semantic modeling strategies. The lexical component captures orthographic patterns and contextual cues characteristic of SW/DB names, while the semantic component leverages contextual embeddings generated by transformer-based language models such as SciBERT, combined with an explicit token-masking strategy to enhance entity-focused representations. A large training corpus was constructed automatically through a hybrid pipeline that integrates citation-hinted extraction with large language model-assisted distillation. Evaluation on two independent benchmark datasets and real-world research articles demonstrates that SNAIL substantially outperforms existing approaches, including domain-specific methods such as bioNerDS2 and general-purpose large language models such as ChatGPT, Gemini, Grok and Claude. Applying SNAIL to large-scale literature analysis further reveals distinct journal-level preferences across bioinformatics subfields. These results demonstrate that SNAIL provides an accurate and scalable solution for identifying bioinformatics resources in scientific texts and enables systematic meta-analysis of tool usage and research trends.

High-phenolic plant genera present a major technical limitation in genomic research. Standard extraction approaches that perform reliably across diverse flora often perform poorly when applied to recalcitrant taxa, producing low DNA yield and integrity incompatible with sequencing requirements. The genus Anigozanthos (Kangaroo paws) from the family Haemodoraceae exemplifies this problem. We identified key physicochemical factors governing extraction failure in this genus and resolved them through targeted modifications to lysis chemistry and contaminant management. The resulting protocol achieved a near threefold improvement in DNA purity, substantially reducing contaminant carry over and consistently yielded high-integrity, long DNA fragments (DIN > 7) across a diverse sample set spanning cultivated and wild material across four diverse genera of Haemodoraceae. We also tested a straightforward purity assessment framework that can be implemented in any standard molecular laboratory, enabling rapid pre-submission quality assessment without the need for specialised equipment. Together these advances open a practical path to genomic characterisation of Anigozanthos that establishes a transferable model for genomic research across Australia ' s chemically complex native flora.

Amplicon sequencing of the 16S rRNA gene is the most widely used approach for profiling bacterial communities, but its taxonomic resolution is typically limited to the genus level. Many species carry multiple divergent 16S rRNA alleles that overlap across species boundaries, an ambiguity that even full-length, long-read sequencing cannot fully resolve. Shotgun metagenomics achieves species-level resolution but remains costly, particularly when only a single genus is of interest. Amplicon sequencing of rapidly evolving, protein-coding housekeeping genes offers a cost-effective alternative, yet no tool exists to identify suitable primer sets for a given target taxon. Here we present AmPair, a Snakemake pipeline that, given a target genus and one or more candidate housekeeping genes, designs and ranks primer pairs binding conserved regions while flanking a variable region capable of species-level discrimination, and validates them in silico across all available genomes. Using the genus Bacillus and the housekeeping gene tuf as a case study, the primer set recommended by AmPair amplified 99% of 2,392 genomes; only 0.04% carried multiple alleles and none showed inter-species allele overlap, compared with 91.41% and 69.49%, respectively, for the standard 16S rRNA V1-V9 region. Applied to a Bacillus community profiled by Nanopore sequencing, the same primers resolved closely related species. AmPair thus offers a generalizable and accessible route to species-level community profiling.

Causal mediation analysis is widely used to identify biological pathways linking exposures to outcomes, but most methods assume homogeneous mediation effects across individuals. In high-dimensional omics settings, this assumption can mask important heterogeneity driven by demographic, genetic, or environmental factors. We propose the M-high-learner, a flexible framework for detecting heterogeneous mediation effects with high-dimensional mediators. The method identifies mediators with subgroup-specific indirect effects while distinguishing them from null or homogeneous signals and controlling the type I error rate. It is computationally efficient, scalable, and yields interpretable sub-types. Simulation studies show that the proposed approach achieves high power while maintaining accurate error control. Applications to the Framingham Heart Study and the Multi-Ethnic Study of Atherosclerosis reveal that the mediation role of gene expression in sexs effect on high-density lipoprotein varies across subgroups defined by body mass index and age. Our framework provides a practical tool for uncovering heterogeneous biological mechanisms in high-dimensional genomic studies. Author SummaryBiological processes linking risk factors to disease often differ across individuals, but many existing methods assume these processes are the same for everyone. This can hide important differences between groups. We developed a powerful method to identify when these pathways vary across subgroups using large-scale molecular data. Our approach detects differences in how intermediate biological factors contribute to outcomes in populations defined by characteristics such as age and body mass index. Applying our method to population studies, we found that some biological pathways operate differently across groups, suggesting that key mechanisms may be missed when differences are ignored. Our work provides a tool to better understand how disease-related processes vary across individuals, which may support more targeted and personalized approaches to health research.

Proteins are dynamic molecules existing in diverse conformational states underlying their biological functions. Although recent approaches have enabled diverse conformational sampling by emulating molecular dynamics simulations, perturbing evolutionary information, or steering internal mechanisms of structure prediction models, predicting conformations resulting from major domain motions or motions that occur over long timescales still remains a challenge. To this end, we introduce GNMCADS, a conformational sampling strategy that enhances the diversity of protein diffusion models by selectively annealing the conditioning signal guided by the intrinsic dynamical organization of the sampled protein. Further, we implement GNMCADS in the diffusion module of AlphaFold3, enabling the generation of diverse protein conformations. When benchmarked across 92 proteins that include 54 class A GPCRs, 15 transporters, and 23 proteins with major domain movements, GNMCADS exhibits improved sampling diversity compared to other current conformational sampling methods.

Network-based analyses of molecular interactions are useful for interpreting high-throughput omics data and identifying therapeutic targets. Cytoscape is the standard platform for these tasks, but users face a trade-off between accessible graphical workflows that are difficult to document and reproducible automation in Python or R that requires programming expertise. General-purpose coding assistants can generate Cytoscape Automation scripts, but remain external to Cytoscape. We present CyChat, a Cytoscape Desktop app that integrates a chat interface and a large language model (LLM) agent into the application. CyChat translates natural language into executable Cytoscape Automation workflows, runs generated Python code, and exports chat sessions with executed code as standalone Jupyter notebooks. To reduce setup barriers, CyChat includes an embedded Python runtime and supports both cloud-based and locally hosted LLMs. CyChat was evaluated across ten Cytoscape workflows using seven LLM providers, each represented by one LLM. The strongest configuration achieves a pass rate above 99%. In a qualitative evaluation based on a published network visualization, CyChat completes the task in 1.5-5 minutes, compared with 15-20 minutes for manual GUI workflows by computational biologists. CyChat is available through the Cytoscape App Store at https://apps.cytoscape.org/apps/cychat.

Bacteriophages can only be understood through their interactions with bacterial hosts. As environmental sequencing efforts expanded, the number of available phage genome sequences has exploded, yet the vast majority of these sequences lack host information. Predicting the host of a newly observed phage is therefore a key challenge in virology. Several computational tools can predict phage-host relationships from genomic data, but they share notable limitations: (1) the number of different hosts that can be predicted remains relatively restricted; (2) tools tend to assign confident host predictions to non-viral input sequences; and (3) most tools have a trade-off between accuracy and speed. Here we present PhageTransformer (PT), a deep learning model for phage-host prediction that addresses these limitations. We benchmark PT against existing tools on 3,881 independent phage-host pairs from GenBank and public HiC data, and demonstrate that it achieves competitive or superior prediction accuracy at greatly reduced runtime.

Scaling the structural annotation of protein-coding genes to all eukaryotic genomes remains a major challenge. While recent deep learning methods rival evidence-based pipelines without requiring RNA-seq or alignments, they are entirely supervised. They depend on large, high-quality training sets from diverse genomes, leaving many basal eukaryotic clades without an accurate ab initio gene finder. We present Vipsania, the first unsupervised deep gene finder. A differentiable hidden Markov layer inside a deep sequence model learns to predict gene structures from unannotated genomes alone. Vipsania is pretrained for virtually all eukaryotes and finetunes without supervision on the target genome. It is, on average, more accurate than supervised methods across most clades and avoids the accuracy drop that supervised models suffer on distant target genomes. Vipsania adapts to non-standard genetic codes and provides a fast and highly versatile tool for unbiased, pan-eukaryotic genome annotation. The source code is available at https://github.com/gaius-augustus/vipsania.

Abstract Imaging-based spatially resolved transcriptomics (imSRT) technologies provide high-throughput molecular-resolution spatial characterization of genes within cells. Conventional analysis methods to identify cell-types and states in imSRT data rely on gene count matrices derived from tallying the number of mRNA molecules detected for each gene per segmented cell, thereby overlooking subcellular heterogeneity that can be useful in defining cell states. To take advantage of the molecular-resolution information in imSRT data and potentially identify cell-states based on subcellular heterogeneity, we developed STARIT (Spatial Transcriptomics As Rasterized Image Tensors). STARIT converts transcripts within segmented cells in imSRT data into an image-based tensor representation that can be combined with deep learning computer vision models for downstream analysis. Using simulated and real imSRT data, we demonstrate that STARIT distinguishes transcriptionally distinct cell-types and further separates cell states based on subcellular transcript localization, which conventional gene count analysis fails to capture. By providing a standardized framework to encode subcellular molecular information in imSRT data, STARIT will enable deeper insights into subcellular heterogeneity and enhance the identification and characterization of cell-types and states that are overlooked by gene count representations.

Recent advances in protein structure prediction, exemplified by AlphaFold, have largely addressed the determination of static structures, one aspect of the protein folding problem. However, predicting folding pathways, by which proteins reach their native states, remains a significant challenge. Here, we present PathFold, a deep learning framework that predicts protein folding pathways directly from sequence information. PathFold leverages an AlphaFold-based module to extract structural information from the sequence and generates a progressive folding trajectory from an extended conformation using a diffusion model. By modeling the full trajectory, it enables prediction of folding intermediates and transition pathways, analogous to those observed in steered molecular dynamics (SMD) simulations. The predicted pathways reveal well-defined intermediates and sequential folding events, and show agreement with experimental folding data, including measured {Phi}-values.

Improvements in data availability, sharing, and integration, together with the development of explainable artificial intelligence (XAI) techniques, are advancing precision medicine for pediatric cancer by facilitating diagnosis, biomarker discovery, and drug development. Data sharing commons and initiatives like the Childhood Cancer Data Initiative (CCDI) provide access to pediatric-specific genomic and clinical data cohorts and improve data availability for pediatric cancer research. Based on CCDI, a scalable AI platform, Graph Artificial Intelligence for Pediatric Oncology (GAIPO), integrates various data modalities from bulk and single-cell omics data to clinical information. Such multi-modal data facilitates the training and development of advanced XAI models for pediatric cancers. We then developed an end-to-end multi-modality framework, PCGS, for pediatric cancer by incorporating omics-specific representation learning via GNN models with cross-attention fusion and multi-objective learning for downstream tasks such as classification, clustering, and survival analysis. This framework outperforms previous supervised multi-omics integration baseline approaches based on glioma and Wilms tumor cohorts and enables GNN model explainability via Shapley value-based feature attribution approaches to explain the contributions of gene-level features across various biomedical tasks, including classification and survival. Given specific background samples (e.g., age groups, sex, grades) as baselines, this explainable GNN model estimates and ranks the importance scores for input features from each omics modality. It identifies background-specific key features for biomarker discovery, risk group identification, and survival analysis in glioma and Wilms tumor, with potential applicability to other pediatric cancers.

Although the biosynthetic pathways of benzylisoquinoline alkaloids (BIAs) have been extensively investigated in several plant species, their transcriptional regulatory mechanisms remain only partially understood. Jasmonate (JA)-responsive group IX APETALA2/Ethylene Responsive Factor (AP2/ERF) transcription factors (TFs) are well-known regulators of specialized plant metabolism, including the biosynthesis of various alkaloids. However, their specific roles in BIA biosynthesis remain largely elusive. Here, we isolated five novel group IX AP2/ERF TFs, designated Benzylisoquinoline alkaloid Jasmonate-responsive AP2/ERF (BJE1-5), from Coptis japonica. Phylogenetic analysis revealed that Benzylisoquinoline alkaloid Jasmonate-responsive AP2/ERF (BJE) proteins belong to subclades distinct from group IXa, which contains well-known AP2/ERF TFs involved in alkaloid biosynthesis. Transient expression analyses in C. japonica protoplasts demonstrated that certain BJEs, particularly CjBJE3 and CjBJE5, positively regulated BIA biosynthetic genes through a mutual regulatory network among BJE members. Moreover, CjBJE3 expression was regulated by CjbHLH1, a unique-type basic helix-loop-helix (bHLH) TF specific to BIA-producing plants. Furthermore, heterologous expression of CjBJE3 and CjBJE5 in cultured Eschscholzia californica cells significantly enhanced the overall BIA production, particularly by increasing end-product benzophenanthridine BIAs, highlighting several uncharacterized biosynthetic genes clustered in the genome. Our findings suggest that BIA-producing species have developed a specific regulatory network comprised of CjbHLH1 and BJE TFs, providing valuable clues for identifying novel biosynthetic enzymes.

MicroRNAs (miRNAs) are key post-transcriptional regulators of development, stress response, and secondary cell wall formation in woody plants, yet annotations for Eucalyptus grandis, the world's most widely planted hardwood, remain fragmented across studies using incompatible discovery pipelines and filtering criteria. Here we present the Eucalyptus MicroRNA Archive (EMA), a curated, locus-resolved database integrating three independent small RNA sequencing datasets spanning vegetative tissue, somatic embryogenesis, and mechanically induced tension wood formation. Applying annotation criteria aligned with current plant miRNA standards, EMA catalogs 99 curated miRNAs (31 previously described, 68 novel) organized into 34 family-level groupings under a three-tier confidence system, known-reference-supported, multi-study replicated, or single-study, that preserves study-of-origin and sample-level evidence for every entry. Cross-study comparison showed that only 9 of 99 entries (9.1%) were independently supported by all three datasets, supporting an evidence-tiered rather than binary annotation scheme. Target prediction against the E. grandis transcriptome yielded 1,773 miRNA-target interactions spanning 764 loci, integrated into a combined miRNA-target and protein-protein interaction network. This network resolved into functionally coherent, mutually isolated clusters, including an miR482-associated NBS-LRR/TIR disease-resistance hub with a substantial translational-repression component, alongside modules enriched for ribosome biogenesis and translation, DNA replication, and nitrogen and carbohydrate metabolism. EMA is publicly accessible through an interactive web dashboard, with all curated data, source code, and analysis scripts openly available, providing a reproducible, extensible framework for E. grandis miRNA research and a template for similarly structured resources in other non-model woody species.

Metabolic dysfunction-associated steatotic liver disease (MASLD), the leading cause of chronic liver pathologies worldwide, represents a growing clinical burden. Its diagnosis remains reliant on liver biopsy that limits early detection and the ability to capture molecular changes across disease progression. A systematic understanding of stage-dependent gene expression changes is essential to identify biomarkers and effectively characterize disease mechanisms. Therefore recent studies provided databases for searching genes as well as prediction of multi-gene signatures for disease progression. However, there is still a need for interactive and comprehensive meta-analysis of datasets of MASLD patients with available histological metadata. Herein, we performed a meta-analysis of RNA-seq datasets using NAFLD Activity Score (NAS; n = 897) and fibrosis stage (n = 856) upon conducting pairwise comparisons across histological stages and identified differentially expressed genes associated with disease progression. Most importantly, we provide our findings via a dedicated web server, the MASLD-META NETWORK (https://masld.scilicium.com), enabling users to interactively explore meta-analysis results across diverse network modalities. In addition, we characterized gene expression dynamics across increasing disease stages to identify consistent progression-associated pathways using Louvain clustering. Network-based parameters such as centrality in combination with meta-analysis scores further highlighted central genes and pathways implicated in disease mechanisms. Accordingly, MASLD-META NETWORK enabled an integrative reassessment of recently published gene signatures, identifying COL1A1, COL3A1, THBS2, FBLN5, and PDGFA as the most central genes, and SULF2, MMP14, IL32, GPNMB, and COL3A1 as candidate markers of earlier transcriptional alterations. Network analysis of MASLD associated biological modules further identified LAMA2 and LAMA3 as previously unrecognized central candidate targets.

Three new species of sun orchid (Thelymitra) endemic to Aotearoa New Zealand are here described. These are T. palustris, T. scabrifolia and T. semaphora. The morphological distinctiveness of these three species has been acknowledged for decades; however, their taxonomic status has remained unresolved. Evidence from existing karyological data, recently generated DNA sequence data (LFY and ycf1) and morphological studies from historical and fresh collections are used here to support their formal description. Both, T. palustris and T. semaphora are restricted to wet habitats north of Auckland (North Island). Thelymitra scabrifolia inhabits mostly scrub, and it has a similar northern North Island distribution, but is has been found also in Manawat[a]whi / Three Kings Islands and historically in Otago (South Island). All three species are polyploids and are of conservation concern.